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1. Introduction, 

Some image analysis tasks depend on the availability of a registered sur- 
face model. Registration can be accomplished using manually identified ground 
control points or by matching the real image with a synthetic image calculated 
from the surface model using assumed reflectance properties. In either case, 
the form of the transformation from image coordinates to model coordinates must 
be known. The registration process is then used to determine the unknown par- 
ameters of the transformation. 

We show here that in the case of satellite images obtained by a mechanical 
scanning system, such as that used on the LANDSAT satellites, an affine trans- 
form applies, if small, second-order effects are neglected. Such a transforma- 
tion has six parameters which depend on the state of the scanning platform. 
Each parameter is exhibited as a function of the components of this state and 
other relevant fixed quantities. These equations can then be used to predict 
transformation parameters if the state of the scanning platform is known. 

A possible application of automatic registration of images and surface 
models is the determination of the parameters of a satellite's orbit. Unfor- 
tunately, a rigid body has six degrees of freedom (position and attitude) and 
so its state has twelve components (position, velocity, attitude and attitude 
rate). Clearly, then, knowing the six parameters of the affine transformation 
at one instant of time is not sufficient to permit a calculation of the vehicle's 
state. 

A series of determinations of the transformation for images taken of dif- 
ferent areas of the earth, however, may permit determination of the vehicle's 
orbit. If we ignore small perturbations, then the position and velocity of 
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the center of mass of the vehicle at one instant of time fully determine its 
orbit. We prefer to use orbital parameters to describe this component of the 
state in some circumstances. 
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2. Basic Orbital Geometry. 

LANDSAT is in a near-polar, retrograde, sun-synchronous orbit which is 
nearly circular. The nominal parameters of this orbit include a semi -major 
axis of 7,294,690 meters, that is, 916,525 meters above an earth with equational 
radius of 6,378,165 meters and oblateness of 1 in 298.3. The nominal period is 
103.267 minutes, which brings the sub-satellite point back to the same spot on 
earth after 251 orbits in 18 days. At the equator, neighboring sub-satellite 
tracks are spaced 159,380 meters. The descending node is nominally passed 
at 9:42 A.M. The orbital inclination is nominally 99.092°, which brings the 
satellite within e = 9.092° of the north pole at the vertex V (see Figure 1). 
All of the parameters drift with time due to perturbing influences such as 
/- v solar wind, light pressure, atmospheric drag, non-spherical distribution of 
masses in the earth, effects of mass expulsion, and so on. The orbit is re- 
adjusted at time using small gas discharges to maintain the positions of the 
ground- tracks within about 35 km of nominal and to prevent the time of north 
to south crossing of the equator from drifting too far from the nominal 9:42 
A.M. The orbital data is derived from radio tracking information. 

Points on the orbit may be conveniently referenced to the vertex V. The 
orbital travel distance p is measured from it, and the reference meridian 
passes through it. The change in geographical longitude X is measured from 
the reference meridian (see Figure 1). 

Ignoring for a moment the rotation of the earth, we find that the nominal 
position of the satellite, S , lies at (geocentric) latitude <j>\ The nominal 
heading of the satellite relative to the local meridian is given by the angle 
H $ . The relationship between orbital parameters e, p and the geographical co- 
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ordinates a $ , <j>' can be established using products of rotation matrices [1]. 
Here we follow a more direct route using spherical trigonometry. 

Considering the right spherical triangle N E g S g (see Figure 1), applying 
the sine theorem, one finds that 



sin (<f>')/sin (90° - e ) = sin (90°- P )/sin (90°) 



or, 



sin <f>' = cos e cos p (1) 

^ Next, from right spherical triangle V P S , one finds, 



sin a /sin p = sin H/sin e 
^ s 



and, applying the cosine theorem (for angles), 



cos A s = - cos H s cos 90° + sin H $ sin 90° cos p 



or s 



cos A = sin H cos p 



Hence, 

tan A s = tan p/sin e (2) 
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Equations (1) and (2) determine geographical coordinates <f>', a in terms of 
orbital parameters e, p. Similarly, 



or, 



Hence, 



cos H_ - - cos A c cos 90° + sin A„ sin 90° cos e 
s s s 



cos H = sin A cos e 



tan H $ = tan e/sin p (3) 



Equation (3) determines nominal heading H s in terms of orbital parameters e, p. 

As it turns out, the earth does rotate and so the sub-satellite point is 
displaced an additional amount in the direction of the local geographical para- 
llel, from point S $ to point S. The latitude $' remains unchanged, of course, 
while the longitude is increased by X Q and the sub-satellite track deviates 
by an angle H e from the nominal direction. In order to calculate these quantities, 
one must know the angular velocities of the earth and of the satellite in its 
orbit. Let these quantities be oa and o> , respectively. 

Since the satellite retraces its path almost exactly every 18 days, after 
completing 251 orbits, we know the ratio of these two quantities, 

O r w = 03 e /a) s = d A e /d P = 18/251 (4) 
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Then, 



A e = r ,/ (5) 
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Actually, w is not constant, unless the satellite is in a circular orbit — 
we will return to this point later. 

At latitude $', the earth surface is displaced a distance Ra> cos <f>' dt 
in a time interval dt, where R is the distance of the surface from the 
earth 1 s center. 

The calculation of the change in heading is a little bit more complicated. 
If we let the satellite heading H = H e + H , then one can see that (Figure 2) 



tan H = (r cos <j>' + sin H $ )/cos H (6) 



Next, 

tan H e = tan (H - H s ) = (tan H - tan H $ )/(l + tan H tan H ) 

tan H = r cos (j) 1 cos H c /(1 + r cos <fr' sin H ) (7) 

"CO o CO S 



Now, from the right spherical triangle PVS (Figure 1), one obtains 



sin (90° - (f,')/sin (90°0) = sin P /sin H g 



That is, 

cos (j) 1 sin H = sin e (8) 
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Similarly, one finds, 



sin (90° - <j>')/sin (90°) = sin P /sin x g 



That is. 



cos <j>' = sin p/sin x (9) 



In deriving (3), we determined that cos H = sin x cos e, so that 



cos <f>' cos H s = sin p cos e (10) 



Finally, we can use equations (8) and (10) to simplify the expression for 
H e (7), 



tan H e = r cos e sin p/(l + r sin e) (n) 



w 



Equation (11) determines H g in terms of orbital parameters e, p and the constant 
r^. Note that r^ sin e is quite small (.0112) and can be ignored in approximate 
calculations [1]. 

Also, now using (8) and (10), we can simplify the expression for tan H (6), 



tan H = (r cos 2 <*>' + sin e)/(cos e sin p) 



or, 



tan H = [rjl - cos 2 e cos 2 p) + sin e]/[cos e sin P ] (12) 
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To sum up, given e and *', we find the orbital travel distance p using 

equation (1), the longitude relative to the reference meridian X = x + x 

s e 

using equations (2) and (5) and the heading H = H g + H using equations (3) 
and (11), or (12). 
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3. Local Solar Time at the Point of Observation. 

An immediate application of the results developed so far is the determina- 
tion of the local solar time at the sub-satellite point. (This gives one some 
idea of what the position of the sun is likely to be). Let the time of the 
descending node be T . That is, the satellite crosses the equator from North 
to South when the local solar time is T Q (for LANDSAT this is nominally 9:42 A.M. 
but tends to vary as the orbit drifts and is readjusted). 

If a point is observed when the satellite has progressed p in its orbit 
from the vertex V, then it still has to travel through an angle (90° - p) be- 
fore reaching the equator. This will take an amount of time which can be ex- 
pressed in hours as r (90° - p)/15. 

Furthermore, the point of observation lies (90° - x ) ahead of the point 
of equator crossing in longitude. Thus the local solar time is later by a time 
which, when expressed in hours, comes to (90° - X )/15. Finally, then, one 
sees that the local solar time at the sub-satellite point is 

T = T Q + (90° - X s )/15 - r u (90° - P )/15 (13) 



where tan x = tan p/sine (2). Because r = 18/251, one finds that the first 
term predominates. As a result, points North of the equator, imaged earlier 
in the orbit , are observed at local solar time after T , while points South 
of the equator, imaged later in the orbits, are observed at local solar time 
before T . 
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4. The Scanning Platform. 

The satellite uses an oscillating mirror to produce the across-the-track 
scan. Individual lines of the image are obtained by this means. The satellite's 
motion in orbit provides for the other scanning direction. Successive lines 
are displaced along the sub-satellite track. Nominally, the optical system 
points straight down and the mirror scanning motion is perpendicular to the 
velocity vector of the vehicle. In practice, there are small but significant 
departures from this ideal state (Figure 3). 

Pitch and roll are measured to an accuracy of .07° using horizon scanners 
sensitive to the infrared radiation (around 14 pm) emitted by the atmosphere. 
Yaw is measured with similar accuracy using a gyro compass. Pitch and roll are 
maintained within ± .4° using the vehicle's attitude control system, while 
yaw is maintained within ± .7°. A major component of the attitude control 
system is a set of inertia wheels which are used in order to keep down gas 
expenditure. 

An attempt is also made to minimize rates of change of attitude which re- 
sult from adjustments. The maximum attitude rates are .015 degree/second. 
Attitude rates are estimated from time-histories of measured attitudes. For 
further information on the scanning platform and its motion, see references 
[1 - 4]. 

Ground tracking information provides good ephemeris data. However, since 
a picture cell in the image is only about 79 meters by 56 meters, one cannot 
expect the position of the satellite to be known accurately enough to predict 
exactly which point of the earth is imaged. Similarly, on-board horizon 
sensors permit a good determination to be made of the attitude of the satellite 
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platform. Nevertheless, these measurements are not accurate enough to permit 
the direct calculation of the ground coordinates corresponding to a particular 
picture cell. Errors of several kilometers may be encountered when this is 
attempted [3]. 

"Precision processing" of satellite image information entails the manual 
identification of known ground control points on each image and the derivation 
of a suitable transformation based on this information. So far, this has 
proved too expensive and LANDSAT images are "bulk processed", that is, treated 
as if the calculated position and attitude of the satellite were exact. As 
a result, the final photographic products may have errors in translation of 
several kilometers. Fortunately, non-linear effects introduced by this ap- 
proximation are small. 

One might envision systems which automatically register image information 
with map or surface model information. In such a system, one has to model 
the imaging operation so that the registration process can be used to deter- 
mine the unknown parameters, such as satellite position and attitude. A clear 
understanding of the scanning process is required to accomplish this. 
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5. Fineness of the Scanner Model. 

A large variety of effects contribute to the imaging transformation. 
Amongst these are large effects which must be considered, such as the motion 
of the satellite in its orbit and the rotation of the earth beneath it. There 
are also smaller effects which have to be judged individually. Some of these 
produce non-linear effects. Examples are panoramic scan distortion (the mirror 
scans evenly in angle, not tangent of the angle), perspective projection (which 
can be dealt with only if a surface model is available) and second-order effects 
of errors in attitude of the spacecraft. The relative importance of these ef- 
fects has already been discussed by others [1 - 4]. The most important criterion 
for including an effect in our model was linearity. 

Fortunately, all major components of the image transformation turned out 
to produce linear transformation of image coordinates. Second order, non-linear 
effects were neglected, but turn out to contribute errors which are typically 
smaller than a picture cell in size. Compounding these linear transformations 
leads to an overall affine transformation which is easy with which to deal. 
Such a transformation has six parameters, which may be found using the registra- 
tion of the image with some surface information in the form of a map or a digital 
terrain model. 

The six parameters, as one might expect, depend rather directly on the 
position of the satellite in its orbit, the attitude of the scanning platform, 
the orbital velocity, and the mirror-sweep velocity. It is conceivable that a 
system which automatically determined the parameters of the affine transformation 
using a matching process of real with synthetic images obtained form a terrain 
model, could also then proceed to estimate the underlying orbital data. A 
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satellite equipped with such a system would be able to determine its position 
or attitude more accurately than it might using predicted ephemeris data ob- 
tained from expensive ground tracking efforts. 
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6. Nominal Parameters of LANDSAT Imaging System (Orbital Parameters Drift), 



Orbit: 



Apogee ~ 917 km Perigee = 898 km 

Inclination = 99.1° (Retrograde orbit) 

Anomalistic period z 103.267 minutes 

[That is, 251 orbits in 18 days] 

Equatorial Earth radius = 6378 km 

Polar Earth radius = 6357 km 

Equatorial speed of rotation = 463.8 m/sec 

Average ground track speed of satellite = 6457 m/sec 



Mirror-Scanner: 



Frequency = 13.260 Hertz 

[That is, 6 lines are scanned every 73.42 msec] 

One line every 12.237 msec 

Lines spaced by = 79.0 meters at nominal height 

390 scans per image 

That is, 2340 scan-lines per image 

[This takes 28.63 seconds and covers = 185 km] 



Pixel Information: 



Instantaneous field of view z 79 m x 79 m 

Mirror amplitude = ±2.886° 

Total scan distance ~ 11.545° 

That is = 185 km at nominal altitude 

Pixels per line (nominal) ~ 3240 

Sampling interval = 9.958 ysec 

That is, about 55.8 - 56.5 m on the ground 

Consequently, u z 6.21 radians/second 

Time to scan (six) lines (in parallel) = 32.238 msec 



Total Image Size: 



2340 x 3240 = 7,581,600 pixels 
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7. Image Coordinate Transformation. 

Let the pixels be numbered sequentially within each scan line and let the 
scan lines be numbered sequentially. Then x $ will be the number of a pixel 
counted from the beginning of a scan-line, while y will be the number of a line 
counted from the beginning of a particular image. (Actually, this is arbitrary 
since the scanner does not start or stop at image boundaries; the continuous 
stream is segmented into images by ground processing). These will be called 
satellite coordinates. 

Now erect a coordinate system in the region of 'interest. First construct 
a tangent plane and let the x-axis run in the west-to-east direction, and the 
y-axis in the south- to-north direction. Now add a z-axis going vertically up 
(we ignore the non-spherical nature of the earth and other such minor effects). 
We will use the notation (x ,y o ) for points on the surface. The satellite 
can also be located in this earth coordinate system. At some reference time 
t Q , it is at (x ,y ,z ) and has attitude a(roll), e(pitch), and y(yaw) (Figure 
3). The three attitude angles will be assumed to be small. 

At time t Q , the scanner will also be at a particular point in its scan 
of the image. Let it be scanning the x -th pixel in the y -th line of the 
image. If the sensor were pointing straight down (that is, a = and g = 0), 
it would be imaging the sub-satellite point (x ,y ) (Figure 4). 

At this point we introduce a convenient artifice, a spherical earth fixed 
relative to the orbit of the satellite. That is, a spherical surface which 
is also sun-synchronous, rotating once a year. Later we will take into ac- 
count the fact that the earth rotates underneath the satellite. We will first 
develop the coordinate transformation for the case of a fixed surface because 
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it is easier to understand this transformation. 

Here it is convenient to refer pixel locations to the reference point 
(x so ,y so ) (Figures 4 and 5). 

x i = (x s - x so } and yi = (y so - y s ) 04) 

Let the angular scanning velocity produced by the mirror during its linear 
phase be u, m (about 6.21 rad/sec) and let t g (9.958 ysec) be the sampling in- 
terval during the scan, then, on a surface at a distance z from the satellite 
and perpendicular to the extension of the optical axis of its scanning system, 
we find a cross-track scanning amplitude x 2 as follows, 

X2 = K Z o V x i (15) 

In the along-the-track direction, the motion of the satellite in its orbit 
provides for the scanning and so, 

y 2 = u s R t £ ) Yl (16) 

where R is the distance of the surface from the center of the earth (about 
6370 km), while u> $ is the angular velocity of the satellite in its orbit (about 
1.014 mi Hi -rad/sec) and t £ is the interval between successive scan-lines (12.237 
mi Hi -seconds). [Actually sU lines are scanned simultaneously eyery 73.42 
mi Hi -seconds.] 

At this point, we note that because of possibly non-zero yaw, the across- 
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track scanning may not be perfectly perpendicular to the along-track scan. 
This skewing effect can be taken care of as follows (Figure 6), 

x 3 = x 2 cos y and y 3 = y 2 . X2 S1 - n y ( 17 ) 

We still have to deal with the effects of roll and pitch. For small angles, 
these will have the effect of shifting the imaged area by an amount proportional 
to the product of the angles and the distance to the surface being imaged. 
Secondary, non-linear effects (such as bending of the scanning line) will be 
ignored, as will non-commutativity of rotations. 

Thus the effects of non-zero roll and pitch can be introduced, 

x k = x 3 - a z q and y h = y 3 - 3 z Q 

where z Q is the height of the satellite above the surface as before. The co- 
ordinate system above lies on the tangent place of the (fixed) sphere. One 
coordinate axis (y) points backward along the sub-orbital track, while the 
other (x) lies at right angles to it. We would prefer to work with a system 
which is aligned with local north. The angle between the local meridian and 
the sub-satellite track (on the fixed earth) is H . We can rotate coordinates 
into a new system as follows (Figure 7), 

x 5 = X4 cos H s + y h sin H g (19) 

y 5 = -xij sin H g + y 4 cos H $ (20) 
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In this new coordinate system, the y-axis points north and the x-axis 
east. Finally, we are ready to introduce the rotation of the earth. It has 
no effect on the value of y, of course, but does introduce a shift in x which 
depends on the time when a particular pixel is imaged. For a particular line 
of the image, this time can be calculated relative to the time t , when the 
reference line was imaged. For line number y s , this time interval equals 
M y s " y so^ = "V 1 ' In tnis t1me interval, the earth has rotated in an 
easterly direction by an amount which depends on the latitude. 

x 6 = *5 + U e R cos $' t £ ) y x and y 6 = y 5 (21) 

^ where a> e is the angular rate of the earth (about 72.722 micro-radians/sec) 
while *' as before is the (geocentric) latitude, and R the distance of the 
surface from the center of the earth. 

To obtain coordinates in the original system (x .y ), we must add the 



e ,J e' 



coordinates of the sub-satellite point (x ,y ) , 

r x oo 



x e = x 6 + x and y e = y 6 + y (22) 
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8. The Overall Transformation. 



All the partial transformations can now be combined, 



0r s 



x = x^ cos H + y k sin H. + (a R t. ) cos <fr' y l + x, 



y e = -x k sin H $ + y k cos H g + y Q 



x e = x 3 cos H g + y 3 sin H $ + (u> e R t £ ) cos <j>' y x - 



(a cos H + 3 sin H) z + x 

s S' 



y P . = ~ x 3 sin H + y 3 cos H - (-a sin H c + 3 cos H) z n + y r 



That is, 



x = cos (H + y) x 2 + sin H_ y 2 + U R t. ) cos <j>' y 1 + x n - 



(a cos H s + 3 sin H ) z 



y = -sin (H + Y ) x 2 + cos H c Y 2 + y n - (-a sin H + 3 cos H e ) z c 



Or, 



x e = cos (H s + y) K z o V* 1 + [sin H s K R V + cos ^'^ e R t s)^ Yl + 



[x Q - (a cos H $ + 3 sin H s ) z Q ] 
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y e = -sin (H s + Y )(w m z Q t g ) Xj + cos H g (a> s R t^) y x + 



[y - (-a sin H g + $ cos H s ) z Q ] 
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9. Form of the Transformation. 

The transformation is of the form, 

x e = a x x + b y 2 + c (23) 

y e = d Xl + e yi + f (24) 

This is an affine transformation, where the six parameters are given by 

a = K z o V cos (H s + y) (25) 

b = U s R t £ ) sin H s + {u q R t £ ) cos <(,' (26) 

c = x Q - (a cos H s + 3 sin H g ) z Q (27) 

d = "K Z o t s ) Sin (H s + y) ( 28 ) 

e = (a) s R t £ ) cos H s (29) 

f = y ~(-a sin H s + 3 cos H $ ) z q (30) 

We can use these equations to predict approximate transformation parameters from 
estimated values of satellite position, attitude and velocity in orbit. Con- 
versely, if we can use ground control points or digital terrain models to 
determine the coefficients of the transformation more precisely, we can try and 
improve the estimates we have of satellite position and attitude. 

From the form of the equations for c and f it becomes immediately clear, 
however, that there are some limits to this process. That is, one cannot dis- 
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tinguish in our model between displacements of the satellite across the track 
and small roll errors. Similarly, displacements along the track have the 
same effects as small pitch errors- Thus two of the six components of position 
and attitude cannot be found this way. 
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10. Attitude Rates, 

Pitch, roll and yaw drift during the scanning of a single image. The 
rates are less than .015 degrees/second. A constant rate of change, 6 of pitch 
has the same effect as a change in along-track ground velocity of 3 z . That 
is, equation (16) becomes, 

y 2 = (co s R + & z Q ) t £ y 2 (31) 

The transformation is altered only in the appearance of U R+ez ) t in 

s o i 

place of (a R tj. Typical values for o> R and hz^ are 6458 and 32 meters/ 

b & 5 

second respectively (when e = .002 degrees/second). This, then, is a small 
f^, but noticeable change (~ .5%). 

A constant rate of change, a, of roll has little effect on the scanning 
of a single line since m m z Q » a z . Successive lines, however, are shifted 
laterally by a z Q t r That is, equation (17) becomes, 

x 3 = x 2 cos y + (« z Q t £ )y! and y 3 = y 2 - x 2 sin y (32) 

This causes additional skewing of the image. 

Here, typical values are w z t = 56.5 meters and a t t = .4 meter 
(for a = .002 degree/second). Again we see a small, but noticeable change 
(~ .7%), resulting in additional skewing (~ .4%). 
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Finally, the transform parameters are now: 

a = (a) m z o V cos (H s + y) (33) 

b = (a) s R + § z Q ) t A sin H g + (a Z() t^) cos H g + (co e R t £ ) cos <(,' (34) 
c = x Q -(a cos H s + g sin H g ) z q (35) 

dss -K z oV s1n ( H s + y) (36) 

e = (a» s R + § z Q ) t £ cos H s - (a z Q t £ ) sin H s (37) 

f = y -(-« sin H s + 3 cos H s ) z Q (38) 

We see here that the transformation parameters depend on the timing (t ,t ) 
of the scanning system, the mirror scan velocity (u> ), the position of the 
satellite relative to the tangent plane coordinate system (x ,y ,z ), the 
orbital velocity vector (as it affects w s and H g ), the attitude of the scanner 
(a, e, y), the attitude rates (a, h, y), and the latitude, <j>'. 

The vertical component of the velocity vector (altitude rate) and also 
the yaw rate do not contribute to linear changes in the imaging transformation. 
The first produces a change of lateral scale from one end of the image to the 
other, the second a tilt of image lines at one end of the image relative to 
those at the other end. Such small, non-linear effects are ignored. Except 
for these two, however, all twelve components of the state of the scanner 
platform influence the transformation parameters. 

If yaw rate and altitude rate are included, one finds small terms in x iyi 
(and y2). The transformation is then no longer an affine transformation. 
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For small regions, the effects of these terms can be ignored -- for images 
which are large fractions of a standard LANDSAT image, they can not. In the 
latter case, one has to include other non-linear terms we have ignored in any 
case and then the transformation can be expressed with sufficient accuracy by 
two second-order polynomials. 
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11. Figure of the Earth. 

To calculate the displacement of points due to the rotation of the earth, 
and to relate the geocentric distance of the satellite to its altitude above 
the surface, one needs to be able to calculate the distance of a point on the 
earth's surface from the earth 1 s center. To a first approximation, a meridonal 
cross-section through the earth is an ellipse (Figure 8) with semi-major axis, 
a = 6,378,165 meters at the equator; and semi-minor axis, b = 6,356,783 meters 
at the poles. 

If we introduce a coordinate system with the x-axis along the semi-major 
axis and the y-axis along the semi -mi nor axis, then the geocentric latitude, 
<£ ! is defined by, 

tan <j>' = y/x (39) 

The more commonly used geographic latitude, <f>, is the angle between a local 
normal to the surface and the equatorial plane. Thus, 

-1/tan * = & (40) 

Using the equation for the ellipse, 

(x/a) 2 + (y/b) 2 = 1 
one finds that, 
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dx K a J y 



so that 



tan ((>' = (b/a) 2 tan <f> (41) 



Further. 



x = (- ab 



/b z + a 2 tan 2 <f,' 



y = ( ab ) tan (j, ' 

^ /b z + a z tan 2 $' 



The distance of a point from the center, then, is, 

^ „. ab 

J /a 2 sirF <j>' + b z cos z <j>' (42) 

The height of a surface feature above mean sea-level must be added to this. 
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12. Orbital Velocity. 

Since one scanning motion depends on the satellite's angular velocity in 
its orbit, it is useful to relate this to other quantities. Using Kepler's 
second law, one immediately sees that 



r dt " 2 (43) 



where r is the radius vector (geocentric altitude of the satellite), while h 
is a constant. Now the area of an ellipse with semi -major axis, a, and semi- 
minor axis, b, is irab, so that 



h 2tt , 
2 = — ab 



where T is the complete period of the satellite. Using w for the angular 
velocity of the satellite, one finds, 

„ _ 2it ab _ . ab , , 

u S -T^- h 7 (44) 

where the average angular rate h = 2ir/T. Further, 

h 2 a 3 = y = GM (45) 

by Kepler's third law, where G is the gravitational constant and M is the mass 
of the earth, y z 396.08 x 10 12 radian 2 meter 3 /second 2 . 
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For orbits with small eccentricity, there is very little difference be- 
tween a and b, so perigee, r . and apogee, r , are given instead, where 

P a 



r p = a (1 - e) 



r g = a (1 + e) 



(46) 



r^. 



and 



Thus, 



a 2 (1 - e 2 ) = b 2 



a = (r a + r p )/2 



b = /r r 
a p 



(47) 



e = 



_§ E. 

a p 



(48) 



Note that perigee and apogee are frequently given as distances from the surface 
of the earth and the equatorial radius of the earth (6,378,165 m) has to be 
added to this in order to find r and r fl . For a more detailed analysis, see 
Lyndanne's modification of Brouwer's analysis of satellite orbits. 
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13. Map Distortion. 

So far we have ignored the fact that the spherical surface of the earth 
cannot be represented on a planar map without distortion. Up to this point, 
coordinates have been referred to a hypothetical plane tangent to the earth's 
surface at a point in the region of interest. Typically, a digital terrain 
model will be derived from a map with a different projection and a transformation 
must be established between the two coordinate systems. 

It can be shown that for typical map projections such as transverse 
Mercator and conformal oblique axis cylindrical there exists a small rotation 
H m of map coordinates, where 



sin H z sin <j> sin (e - e ) 
m o K o y 



(49) 



Here the projection is centered on a point at longitude e and latitude <f> 
and the point of interest is at longitude and latitude $. Consequently, 
the map coordinates, (x m »y m )» are related to the geographical coordinates on 
the tangent plane (x ,y ) by 



m 



r m 



cos H 



m 



-sin H 



m 



sin H, 



m 



cos H 



m 



(50) 



r^ 



where s is the scale of the given map. There will also be a small scale change 
which varies with 1/cos (e - ej for transverse Mercator and with 1/cos U - & ) 

< \ T T Ql 

for conformal oblique axis cylindrical projection. Typically, this effect is 
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so small that it may be ignored. 

The affine transformation, (23) and (24) must be pre-multi plied by an 
augmented rotation matrix M to correctly relate satellite image coordinates 
to map coordinates. 



M = - 
s 



cos H -sin H 

m m 



sin H cos H 
m m 



(51) 



r^ 



/->. 
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14. Example of Map Rotation. 

The national map of Switzerland is based on a conformal oblique axis 
cylindrical projection with center at the city of Bern. 



r^ 



e Q = 7°26'20" 



4> = 46°57 , 10" 



The region of interest covered by an available terrain model lies at 



6 ~ 7°8" 



* = 46°15' 



f*\ 



The map rotation then is -.225° (-.00393 radians). The transformation matrix 
becomes 



1.0 +.00393 
-.00393 1.0 



(The scale error is less than one part in 10,000 and can be ignored). For more 
details regarding suitable map transformations see Colvocoresses paper on the 
"Space Oblique Mercator" projection. 
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15. Numerical Example -- Estimating Transform Parameters. 

LANDSAT image number 1078-09555, produced 1972/0ctober/9 shows a region 
of Switzerland including a mountain range called "Dent de Morel es". The image 
annotation data suggest that at the nadir the geographic latitude was 45.9197° 
and the heading 193.11324°. Using equation (41), we see that the geocentric 
latitude of the nadir point is 45.7274° and so, using equation (8), it appears 
that the orbital inclination must have been e ~ 9.11°. 

The altitude is given as 915,724 meters near the region of interest, which 
lies at an average of 1700 meters above average sea level. So z : 914 km. The 
angular velocity of the satellite is slightly above its average rate, 



p OCT 

w s = 24 x 60 x 60 x Ts~ x 1 - 00967 radians/second 



The region of interest lies at geographic latitude <j> = 46.25° (and thus at 
geocentric latitude $' = 46.06°), while the scanner at that time is above 
geographic latitude <j> = 46.40° (that is, geocentric latitude <j>' = 46.21°). 

Further, 

w e = 24 x 60 x 60 radians/second 



cj = 6.21 radians/second 



t s = 91958 usecond 



*£ = 13.62 x 6 Seconds 
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The image annotation also gives, 



a (roll) = -.20370° 



a = -.00160°/second 



3 (pitch) = ,06688 c 



,00109°/second 



Y (yaw) = .23387' 



Y =.00189°/second 



The earth radius (equation 42) near the region of interest is 6,367,081 meters, 
Adding the average elevation above sea level, we get 

R = 6,368,800 meters 



Using equation (1), one finds that P = 43.021° and by equation (8), that H = 
13.226°. Further, (H g + Y ) = 13.460°. Then also, U|n z Q t g = 56.521 meters, 
(w s R + i z Q )t £ = 79.582 meters, 'a z Q t £ = -.312 meter, and m Q R t = 5.667 
meters. 

So, finally, 



a = 54.969 
d = -13.156 



b = 21.837 
e = 77.543 



c = x Q + 2919 



f = y Q - 1782 



See Figure 9 for a graphical illustration of this transformation. It shows as 
a parallelogram the region of the surface scanned when an image with a number 
of lines equal to the number of pixels per line is gathered. 

We can take the inverse of the 2 x 2 sub-matrix and obtain, 



where, 
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x 1 = .01704 x„ - .00480 y + x' 
e e o 



y x = .00289 x + .01208 y + y' 
1 e ■'e •'o 



y" 



a b n 
d e 



-1 



c i 



/^""s 



See Figure 10 for a graphical illustration of this transformation. It 

shows as a parallelogram the appearance in the image of a square region on the 

ground aligned with the north-south and east-west axes. Finally, if we have 

a terrain model on a 100 x 100 meter grid, such that x = 100 * i and y = 
100 * j, then, 



x x = 1.704 i - .480 j + x' 
yi = .289 i - 1.208 j + y' 

Finally, we have to introduce the map distortion by post-multiplying by the 
inverse of the map transformation matrix introduced earlier (the inverse of a 
rotation matrix equals its transpose). The matrix then becomes 



/^\ 



1 . 702 



,294 



-.487 



1.207 
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16. Determining Orbit Parameters from Transforma t i on . 



Many parameters appear in the equations for the six coefficients of the 
affine transformation. Some are known accurately, others only approximately. 
For example, t g and t^ are fixed fairly accurately by electronic oscillations 
on board the satellite. The angular rate of the earth u is constant and ac- 
curately known, while u> s , the angular rate of the satellite, depends on its; 
altitude, semi-major axis, orbit eccentricity, and orbital period. The angular 
rate of the scanning mirror, u , varies somewhat during each scan, though in 
a fairly repeatable fashion. The earth radius, R, can be calculated with 
sufficient accuracy if the latitude, <f>', is known approximately. 

It is reasonable, then, to assume that one can accurately predict a value 
^ for u e R t cos <f>'. Then let 



b' = b - a) R t cos <f>' 



The following equations permit the determination of useful satellite parameters 
from transformation parameters: 

H s + Y = tan _1 (-d/a) (52) 

K z o V = ^ 2 ~ T ^ (53) 

• la) s R + B z )* + (a z )2t = >/(b'F + e* (54) 



<-*s If the heading H g can be estimated, one finds 
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(« s R + e z Q ) t £ = b 1 sin H s + e cos H g (55) 

(& z ) t £ = b' cos H g sin H g (56) 

Alternatively, the heading can be calculated if the attitude rates are known, 

b' ( u R + e z ) - e (aZ ) 
H s = tan b' (d z p ) f e(" s R -M z°) (57) 

Clearly there are limitations to this, since the full state of the scanning 
platform cannot be ascertained from six parameters alone. A series of such 
measurements is needed to determine all twelve components of state. 
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17. Numerical Example. 

Suppose the transformation matrix applicable to a 100 x 100 meter grid 
was found by image registration techniques using synthetic images, 



1.694 



.300 



.512! 



1 .216 J 



for any area with a map distortion as defined previously. Post-multiplying 
by the map transformation matrix gives 



/-\ 



1.693 



,295 



-.505 



1.217 



The area is at geographical latitude <j> = 46.25°. The geocentric latitude is 
then $' =46.06° and so R = 6,367,081 meters. If the region of interest lies 
at an average altitude of 1700 meters, and t £ = 1/(6 x 13.62) seconds, then 

u e R t % cos <j>' = 3.932 meters/scan-line 

The inverse transformation matrix then is 



55.08 



22.86 



•13.35 



76.63 
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From this one finds, 



H $ + y = 13.624 c 



w m z o *s = 56 - 6 ^ meters/pixel 
And, if the attitude rates are assumed to be very small, 



H $ = 13.874° 



and so. 



^ Y ~ .25V 



and 



U s Rt^) = 78.93 meters/scan-line. 
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18. Using Ground Control Points. 

An alternate method of estimating transformation parameters is based on 
the identification of points of known ground position in the image. Since 
the affine transformation has six parameters, one needs to locate three such 
points. Let the coordinates of the points be (x 1( yj, (x 2 , y 2 ) and (x 3 , y 3 ) 
in the image and (x{, y{), (x 2 , y 2 ) and (x 3 , y 3 ) on the map, then 



xj = a x i + b y i + c 
y\ = d x i + e y i + f 



and so on. Thus. 



Xl 


yi 


i 


x 2 


y 2 


i 


x 3 


ys 


i 



a 




x{ 


b 


= 


x 2 


c 




x^ 



So 



Xl 


yi 


x 2 


Y2 


x 3 


y 3 
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i 
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Similarly, 
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If more than three points can be identified, better accuracy is available 
using a least-squares procedure. That is, if 



M = 



xi yi 
*2 y> 



x n y* ] 



Then a good set of values for the parameters is 



= (M T M) _1 M T 



xi 



X2 



= (M T M) _1 M T 



Yi 

y 2 



Typically, the accuracy obtained by fitting a discrete set of ground control 
points has been found to be inferior to the area-based matching of real and 
synthetic images. 
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